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Using discrete simulations, we investigate the behavior of a model granular material within an 
annular shear cell. Specifically, two-dimensional assemblies of disks are placed between two circular 
walls, the inner one rotating with prescribed angular velocity, while the outer one may expand or 
shrink and maintains a constant radial pressure. Focusing on steady state fiows, we delineate in 
parameter space the range of applicability of the recently introduced constitutive laws for sheared 
granular materials (based on the inertial number). We discuss the two origins of the stronger strain 
rates observed near the inner boundary, the vicinity of the wall and the heteregeneous stress field 
in a Couette cell. Above a certain velocity, an inertial region develops near the inner wall, to which 
the known constitutive laws apply, with suitable corrections due to wall slip, for small enough stress 
gradients. Away from the inner wall, slow, apparently unbounded creep takes place in the nominally 
solid material, although its density and shear to normal stress ratio are on the jammed side of the 
critical values. In addition to rheological characterizations, our simulations provide microscopic 
information on the contact network and velocity fluctuations that is potentially useful to assess 
theoretical approaches. 

PACS numbers: 45.70.Mg, 81.05.Rm, 83.10-y, 8,3.80.Fg 



I. INTRODUCTION 

Significant progress in the modeling of dense granular 
flow in the inertial regime has been brought about by the 
recently introduced viscoplastic laws S|, as identi- 
fied in experiments and discrete numerical simulations 
in two-dimensional (2D) [1, H, 01 and three-dimensional 
(3D) ^MM situations. 

One typically considers homogeneous assemblies of 
grains of size d and mass density pp, under shear stress 
a and average pressure P. Denoting the shear rate as 
7, constitutive laws are conveniently expressed as rela- 
tions between dimensionless quantities: effective friction 
/i* {— <t/P), solid fraction i^, and most noticeably iner- 
tial number / = -yd-^/ Pp/P, thus rescaling various exper- 
imental data into a consistent picture. As the ratio of 
the inertial to shear times, the latter parameter quanti- 
fies the inertial effects. For a frictional material, a small 
value of / (< 10^^) corresponds to the quasistatic crit- 
ical state regime, while a large value of / (> 10^^) cor- 
responds to the collisional regime As / increases, 
solid fraction ly decreases approximately linearly starting 
from a maximum value Vmax = {dynamic dilatancy 
law), while the effective friction coefficient p* increases 
approximately linearly starting from a minimum value 
Pmin — tan0 (dynamic friction law). This yields a vis- 
coplastic constitutive law, with a Coulomb frictional term 
and a Bagnold viscous term. 

In the quasistatic regime (/ 0) this approach in- 
dicates p* — > Pmin constant, independently on the 
strain, in steady shear fiow. Below this minimum 
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stress ratio, quasistatic strains are possible that are de- 
scribed by elastoplastic models. For large enough shear 
strains [ill, a solidlike material approaches the so- 
called critical state, which coincides with the state of 
steady shear fiow in the limit of / ^ 0. 

Once constitutive laws are obtained on dealing with 
homogeneous systems, they should be locally applica- 
ble to all possible fiow geometries. Of course, they are 
quite unlikely to provide a proper description of some 
strongly heterogeneous situations occurring when strain 
is localized near boundaries, in thin layers, on a scale 
of a few grains. Yet, for smoothly varying stress fields, 
they might prove sucessful, as was shown e.g., with fiows 
down inclined planes. Those were studied in the absence 
of lateral walls both experimentally and through discrete 
simulations (see [2| for a review). Then the stress distri- 
bution becomes heterogeneous but the effective friction 
remains constant, so that the situation is comparable to 
homogeneous flows. More remarkably, a three dimen- 
sional version of the constitutive law p, was found to 
model similar flows between lateral walls, which induce 
truly three-dimensional stress distributions and velocity 
profiles [ri |. 

Other simple geometries are the vertical chute and the 
annular shear The present paper investigates the ma- 
terial behavior in the annular (Couette) shear geometry, 
for which the sample is confined between two rough cylin- 
ders and sheared by the rotation of the inner one. The 
annular shear cell is a classical experimental device to 
measure the rheological properties of complex fiuids, and 
has been used for granular materials, both in two dimen- 
sions_fl|Jl6l, It', 18, and in three dimensions [13, [HI, 
[23, [2ll2l^5, 26, 27, 28, 2^ [13, [3l|, [32, 33, sl- In this 
geometry, the stress distribution is well known, as will 
be detailed in the following: the normal stress is approx- 
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imately constant while the shear stress strongly decreases 
away from the inner wall. The decrease of /i* away from 
the inner wall then explains the localization of the shear. 
We may even expect a transition between an inertial flow 
near the inner wall, where fj.* > Umini cind a quasistatic 
regime further, which analysis would help understand- 
ing shear localization near a wall (influence of shear rate 
and conflning pressure on the width and dilation of shear 
bands), of interest in industrial conducts [s^, geotechni- 
cal situations [IBl and tectonophysics js^]. 

Following previous discrete simulations (ssl. Isol. liol. [iH . 
we investigate the rheology and the microstructure of 
granular materials in this geometry. We consider two- 
dimensional, slightly polydisperse assemblies of cohesion- 
less frictional disks. This allows to vary the shear state 
and provides access to microscopic information at the 
scale of the grains and of the contact network, hardly 
measurable experimentally. We prescribe the shear rate 
and the pressure, allowing global dilation of the shear 
cell. To save computation time, we implement periodic 
boundary conditions. All along this paper, we shall com- 
pare our results with the homogeneous shear case [1]. 

Sec. ini is devoted to the description of the simulated 
system, its preparation and the deflnition of dimension- 
less control parameters. Sec. IIIII describes the influence 
of the shear velocity and of the system scale on the shear 
localization near the inner wall through the radial profiles 
of various quantities. Sec. [IV] shows the validity of the 
previous constitutive law for inertial regime, and analyzes 
its limit in quasistatic regime. Sec. [V] then explains how 
the constitutive law is able to predict various quantities 
measured in Sec. IIIII 

Preliminary and complementary results are presented 




II. SIMULATED SYSTEMS 
A. Annular shear 

The simulated systems are two dimensional (2D) as 
indicated in Fig. [TJ The granular material is a dense as- 
sembly of n dissipative disks of average diameter d and 
average mass m. A small polydispersity of ±20% pre- 
vents crystallization. 

The granular material is subjected to annular shear 
between two circular rough walls. The outer wall (ra- 
dius Ro) does not rotate, while the inner wall (radius 
Ri) moves at the prescribed rotation rate VL. The wall 
roughness, which reduces sliding, is made of contiguous 
glued grains with the same characteristics as the flowing 
grains (polydispersity and mechanical properties). We 
call r and 9 the radial and orthoradial directions, r = Ri 
and r = Ro respectively correspond to the centers of the 
grains which compose the inner and outer walls . As the 
grains of the inner wall form one rigid body, the motion 
of each of them combines a translation of its center with 
velocity ftRiCe and a rotation rate fl. 




FIG. 1: Annular shear geometry (black grains constitute the 
rough walls). 



Since we have not observed any influence of the outer 
wall on the behavior of the sheared material close to the 
inner wall for Ro > 2Ri, we have set Ro = 2i?i, in the 
results presented in this paper. The geometry is then 
deflned by the sole value of Ri/d. 

An important feature of our discrete simulations is the 
control of the normal stress exerted by the outer wall on 
the grains, as done in [4l[. We prescribe <7rr{Ro) = P, 
through the radial motion of the outer wall, given by: 
Ro = {P — crrr{Ro))/gp, where gp is a viscous damping 
parameter. In steady state, the motion of the outer wall 
oscillates around a mean value {{Ro) = 0) corresponding 
to a prescribed value of the normal stress at this point 
{{arr{Ro)) = P)- Such control of the radial stress is 
applied in cylinder shear apparatus aimed at studying 
the behavior of soils near an interface [i^. It differs 
from most experiments and discrete simulations, where 
the volume is fixed in two dimensions JS, 
131, or dilatancy ispossible through the free surface in 
three dimensions (Ifl, [13, [H, [H, m, HI, [13, [H, [H . 

We use the standard spring- dashpot contact law de- 
scribed in Q, which introduces the coeflicients of resti- 
tution e and friction ^, and the elastic stiffness pa- 
rameters kn and kt- Discrete simulations are car- 
ried out with standard molecular dynamics method, as 
in 0, H, m, m, [13, [Hi. The equations of motion are 
discretized using Gear's order three predictor-corrector 
algorithm [49| . 

To decrease the computation time, we have introduced 
periodic boundary condition along 9, exploiting the an- 
gular invariance [50|. This reduces the representation of 
the annular shear cell to an angular sector Q < 9 < Q 
{9 < tt) instead of the whole system < 9 < 27r. 
Q = 2Tr/N, where N is an integer. The description of 
this method, together with the analysis used to choose 
the values of 6 according to the size of the systems, is 
presented in App. A. The list of simulated geometries is 
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given in Tab. I. We notice that the studied systems are 
much larger than in previous discrete simulations. 



TABLE I: List of simulated geometries. 





n 


Ri/d Q/2-K 


R25 


1500 


25 1/4 


R50 


3100 


50 1/8 


-Rloo 


8000 


100 1/12 


-R200 


15700 


200 1/24 



B. Dimensional analysis 

In our discrete simulations, the system is completely 
described by a list of independent parameters associated 
to the grains and to the shear state. As a way to reduce 
the number of parameters, it is convenient to use dimen- 
sional analysis, which guarantees that all the results can 
be expressed as relations between dimensionless quanti- 
ties. 

The grains are described by their size d and mass to, 
their coefficients of restitution e and friction ^, and their 
elastic stiffness parameters A:„ and kt. It was shown in 
previous discrete simulations 0, H, [5l|, that kt/kn and 
e have nearly no influence on dense granular flows. Con- 
sequently, kt/kn was flxed to 0.5 and e was flxed to 0.1. 
The influence of /i, especially near ^ = 0, has been shown 
in In this paper, we restrict our analysis to the value 
/i = 0.4, except for the discussion of the constitutive law 
where the case of frictionless grains (/i = 0) will be also 
analyzed. Results for other values of /i may be found 
in |4a |. 

The shear state is described by the prescribed normal 
stress on the outer wall P, the rotation rate of the inner 
wall fi, the radius Ri and Ro of the two walls, and the 
viscous damping parameter gp. We have not observed 
any influence of gp, once the shear zone is localized near 
the inner wall and separated by a relatively thick layer of 
material from the outer wall. The dimensionless number 
gp/^/mkn remains of order 0.1 in all our simulations, so 
that the time scale of the fluctuations of Ro is imposed 
by the material rather than the wall, and that the wall 
sticks to the material. Consequently, the shear state is 
described by the geometric parameters Ri/d {Ro = 2Ri), 
and by the dimensionless tangential velocity of the inner 
wall (also called shear velocity): 




(1) 



which is similar to the notion of inertial number, but 
at the scale of the whole system. A small value of Vg 
corresponds to the quasistatic regime, while a large value 
corresponds to the collisional regime. Seven values of Vg 



have been studied systematically for all systems: 0.0025, 
0.025, 0.25, 0.5, 1.0, 1.5 and 2.5. The value 0.00025 was 
also considered in a few cases. 

Moreover the stress scales kn and P may be compared 
through the dimensionless number k = kn/P. Let us 
call h the normal deflection of the contact (or apparent 
interpenetration of undeformed disks). Being inversely 
proportional to the relative deflection h/d of the con- 
tacts for a conflning stress P, k is called contact stiffness 
number A large value corresponds to rigid grains, 
while a small value corresponds to soft grains. It was 
shown Q that it has no influence on the results once it 
exceeds 10*, which is the value chosen in all our discrete 
simulations (rigid grain limit). 

In the following (both text and figures), the length, 
m ass, t ime and stress are made dimensionless by d, to, 
y^rn/P and P, respectively. 

Table II gives the list of material parameters. 

TABLE IL List of material parameters. 



polydispersity ^ e kt/kn ft 



±20% 



0.4 0.1 0.5 10* 



C. Steady shear states 

For a given sample, the first step consists in depositing 
the grains without contact and at rest between the two 
distant walls. Applying a normal stress at the outer wall, 
we compress the assembly of grains, considering first that 
they are frictionless (^ = 0) , so as to get a very dense ini- 
tial state. Except near the walls, its solid fraction is close 
to 0.85, near the random close packing of slightly poly- 
dispersed disks [H^l- When the granular material sup- 
ports completely the applied normal stress, the grains 
are at rest and the dense system is ready to be sheared. 
We start to shear the material (now considering that the 
grains are frictional) imposing the rotation of the inner 
disk. After a transient, the system reaches a steady state, 
characterized by constant time-averaged profiles of solid 
fraction, velocity and stress. In practice, the stabiliza- 
tion of the profiles depends on the considered variable. 
If we take the inner wall displacement VgAt (where At 
is the simulation time) as a shear length parameter, the 
stresses usually present a short transient on a distance 
around VgAt < 5. However, the stabiHzation of the soHd 
fraction rather requires VgAt sa 50, mostly because of the 
very dense initial state. Consequently we consider that 
the condition to reach a steady state is VgAt > 100. This 
procedure provides an initial state with a shear velocity 
Vg. As a way to guarantee an initial state consistent for 
the comparisons between discrete simulations with dif- 
ferent Vg, the procedure is first applied with the highest 
value of Vg, and then Vg is progressively decreased. 

In steady state, we consider that the statistical distri- 
bution of the quantities of interest (structure, velocities. 
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forces. . . ) are independent of t and 0, so that we average 
both in space (along 0) and in time (considering 200 time 
steps distributed over the distance VeAt > 200). Then 
we calculate the profiles of solid fraction, velocity and 
stress components according to the averaging procedure 
described in App. B. 

Beyond the number of acquisition points, the consis- 
tency of the averaged values depends on the shear strain 
accumulated during the acquisition of data. We concen- 
trate our interest on the region where the system may be 
considered in a steady state, which occurs at large enough 
shear strain. Based on the observation of the transients, 
we consider that this is true when jAt > 10. Because 
of the strain locaHzation, the region of interest is located 
near the inner wall and limited to Ri ^ Ri + Rsteady 
where the value of Rsteady is given in Tab. III. 



Fig. [2^ shows the coarse-grained profiles of the nor- 
mal stress component arr in geometry R^q for different 
wall velocities Ve , while Fig. [2b shows the ratio between 
the orthonormal and the normal stresses — . The nor- 
mal stress (Trr is nearly constant and equal to the con- 
fining pressure P. The darr/dr term in the momentum 
equation ^ smoothes the (!„ profile, which might ex- 
plain the absence of fiuctuations of Grr- A crude es- 
timate of the centrifugal effects may be given, if the 
last term of equation Q is neglected, and, anticipat- 
ing on Sec. IIII Bl and Sec. IIII Cl a constant solid fraction 
V ~ 0.8 is assumed and an exponential velocity profile 
Veil") = Ve exp (— (r — Ri)/£), with £ between 2 and 6: 



~{Ri 



1 |< 



(4) 



TABLE III: Limit of the steady state region. Minimum ano 
maximum values correspond to global quasistatic regime ano 
to Vg = 2.5 respectively. 



-^steady 




7-17 


R50 


9-25 


Rioo 


13-35 


-R200 


18 - 52 



III. LOCALIZED SHEAR STATES 

In this section, we show the shear localization near the 
inner wall through the radial profiles of different quanti- 
ties. In App. C we focus on internal variables associated 
to the contact network [HI, (coordination number Z 
and mobilization of friction M) and to the fiuctuations 
of the motion of the grains, translational or rotational. 
We systematically discuss the infiuence of Ve . 



Stress field 



In steady (^ = 0) annular shear fiows (^ = 0), 
without radial fiow {vr = 0), continuum mechanics pre- 
dicts (H^l a variation of normal stress related to the 
velocity profile, and a decrease of the shear stress 
are associated to the conservation of the torque: 



4i^ Vg 
TT r 



dr 



o'ee 



^(-)^ 

r 



(2) 



(3) 



where iy{r) and vg{r) are the soHd fraction and orthora- 
dial velocity profiles, S is the shear stress at the inner 
wall {S = are{Ri)) and an are positive for compression. 



401), and was previously 
This very small normal 



Consequently, for Ri = 50, | arr{Ri) — 1 |< 0.05 for Ve = 
1 and £ = 5. For Ve = 2.5, the centrifugal effects might 
become significant, however it has not been observed. 

The radial a^ and orthoradial age stresses are nearly 
equal for r — Ri < 10. This has already been observed 
in other configurations (plane shear [1, d, [§] within less 
than 5%, inclined plane [a, |3 
reported in annular shear |40l. l44 1 
stress difference is not explained yet. The fiuctuations of 
aee for r — Ri > lb probably refiect the frozen disorder 
beyond the steady zone, where the material is much less 
deformed than closer to the inner wall, so that the time 
averaging is unsufficient. Consequently these fiuctuations 
increase as Ve decreases. 

The shear stress profiles are{r) shown in Fig. [3K (for 
different shear velocity Ve) are consistent with the 1/r^ 
decrease of The oscillations about the mean value 
are due to the material structuration near the inner wall 
(see Sec. IIIICP and to the frozen disorder in the very 
slowly sheared regions, which is beyond the steady zone. 

Fig.jSb shows the dependence of the shear stress at the 
inner wall S on shear velocity Ve. Below a certain value 
(Ve < 0.025), S tends to a finite limit. Consequently, the 
shear stress profiles cr,-e(r) become independent of Vg. 
This behavior characterizes the global (that is to say, in 
the whole system) quasistatic regime, where the stresses 
(and other state variables) do not depend on the velocity. 
However, for Ve > 0.025, inertial effects become signifi- 
cant and S increases with Ve. Previous works reported 
a similar dependence of the shear stress on the shear ve- 
locity in other configurations (see j2| for a review) . More 
specifically, the experimental measurement of the torque 
as a function of the rotation rate in the annular shear ge- 
ometry indicates a transition from a rate independent to 
a rate dependent regime [13, [13, IIS, S [Hi| . Our results 



can be approximated by a function like S = S, 



qs 



where Sqs is the global quasistatic limit value, a and /3 
are two constants. We notice that f3 is close to 1/2 rather 
than 2 as might be naively expected from Bagnold's rhe- 
ology. We notice that in the experiment of [20| , the tran- 
sition occurs for Ve — 0.3 (after appropriate rescaling). 
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FIG. 2: (Color online) (a) Normal stress crrr{r) and (b) ratio 
between the normal and the orthonormal stresses agg / arr{r) 
profiles for different shear velocities. (W) Ve = 0.0025, (•) 
Ve = 0.025, (m)Ve^ 0.25, (>) Ve = 0.5, (^) Ve = 1.0, (a) 
Ve = 1.5, (k) Ve = 2.5. Geometry Rso- 
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FIG. 3: (Color online) (a) Shear stress profiles t7re(r) (solid 
lines ) and fit according to Eq. ^ ( dashed lines ) for different 
wall velocities Ve- (b) Shear stress at the inner wall S as func- 
tion of Ve (semi-logarithmic scale). The solid line represents 
the function: S = 0.26 + O.lSVg''-^^ Geometry R50. 



which is not far from what is observed in Fig. [3b • How- 
ever, we also point out that the S{Vg) curve, here shown 
for geometry R50, in fact depends on the geometry. 



B. Velocity field 

The shear localization near the inner wall is revealed 
by the strong decrease of velocity profiles vg (r) shown on 
Fig. m The decay appears to be nicely approximated by 
a Gaussian function vg/Vg = exp[~a{r — Ri) — b{r — Ri)'^], 
as shown on Fig. [H We notice however that there is a 
sliding velocity for the higher value of Vg (2.5), which is 
apparent in Fi g. [Sk - Previous studies in 2E) systems (2ll . 
m, [3^, SO, lil. Is?! found an exponential shape, while a 
gaussian decay was observed in three dimensional (3D ) 
systems for non spherical or polydispersed grains |22l |. 
The agreement between the measurement of the velocity 
profiles in 3D experiment (using 3D MRI velocimetry in 
the bulk or CIV at the free surface) [3| and 2D discrete 
simulations is satisfactory [33 |. 



The normalization of vg (r) by shear velocity Vg allows 
to clearly visualize the infiuence of this latter parameter 
on the velocity profiles. In the global quasistatic regime 
(Vg < 0.025), there is no infiuence, while for increasing Vg 
above 0.025, an increase of the localization width is ob- 
served, consistently with experimental observations (25| . 

The shear rate is equal to 7(r) = We 
denote w(r) the profile of the average angular velocity of 
the grains. As previously reported in discrete simulations 
of granular fiows Q, the average angular velocity is equal 
to half the local shear rate (or vorticity) uj{r) = — 7(r)/2. 
Oscillations of the average angular velocity are observed 
in the 3 or 4 first grain layers near the inner wall (Fig.[5|), 
as previously noticed by They may be due to the 

frustration of the rotation of the fiowing grains in con- 
tact with the glued grains of the walls (which rotate with 
angular velocity fl). 
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FIG. 4: (Color online) Influence of the shear velocity Ve on 
the velocity profiles vg(r) (semi-logarithmic scale). (J)Vg = 
0.0025, (U) Ve ^ 0.25, (t^) Ve = 0.5, (^) Ve = 1.0, (<) 
Ve = 1.5, (a) Ve = 2.5. The solid line indicates the function 
ve/Ve = exp[-0.34(r - R,) - 0.0015(r- - 7?,)^], and the dashed 
one the function ve/Vg = exp[-0.21(r- - R^) - 0.002(r - R^f] . 
Geometry Rso- 



C. Solid fraction 

Fig. [6] shows the solid fraction profiles ^(r) for Vg — 
0.025 and 2.5. In the global quasistatic regime {Ve < 
0.025), the profile becomes independent of Vg, while a 
decrease of the solid fraction is observed for increasing 
Vg. The material is significantly dilated near the inner 
wall [l3,[3l,[i3|, and is structured in about 5 layers close 
to the inner wall, with a higher amplitude for low Vg. 
This was previousl y o bserved in various shear geome- 
tries (2^. [sil. [ssl. [sol. l60| . This structuration of the granu- 
lar material certainly affects the sliding of layers of grains, 
with significant consequences on the mechanical behavior 
near the wall. As previously reported [iol. [43|. indepen- 
dently of the infiuence of Vg, solid fraction u increases 
toward a value i^max (close to 0.82, the solid fraction in 
the critical state for frictional disks with a similar poly- 
dispersity away from the inner wall, and remains 

close to its larger initial value 0.85 in the region where 
the material has not been sheared enough. 
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FIG. 5: (Color online) Influence of the inner wall on the col- 
lapse between the average angular velocity iij{r) (hollow sym- 
bols) and the shear rate 7(r) (full symbols) (a) Ve = 2.5 and 
(b)Ve = 0.025. Geometry R50. 



FIG. 6: (Color online) Influence of shear velocity Ve on the 
structuration near the inner wall. Solid fraction profiles v^r) 
(a) in the whole system and (b) in the region close to the inner 
wall. The solid line is an average over 3d, while the dotted 
line is an average over 0.5d. Geometry R50. 
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IV. CONSTITUTIVE RELATIONS 

In Sec, mil and App. C, it has been shown that shear 
velocity Vg, if small enough, no longer influences the 
radial profiles of various quantities (see Fig. [Sj HI [21] 
andES]). Then, the whole system is in the quasistatic 
regime. When Ve increases, the shear rate 7 increases in 
the whole sample. Above a certain level of shear rate, in- 
ertial effects have significant effect on the material behav- 
ior, which characterizes the inertial regime. Considering 
the shear stress distribution in the annular geometry and 
the decay of the velocity away from the inner wall, we 
expect that the inertial zone begins at the inner wall and 
that its thickness increases when Ve increases (Fig. [7]) . 




FIG. 7: Inertial and quasistatic zones. 



In this section, we analyze the relations between dif- 
ferent dimensionless quantities in the inertial regime and 
how they are affected by the transition to the quasistatic 
regime. We restrict our analysis to large enough shear ve- 
locity {Vg > 0.025), so that a wide enough inertial zone 
exists close to the inner wall. 



A. Inertial number and mechanical behavior 

Discrete simulations of homogeneous plane shear flows 
[1| have revealed that the constitutive law of dense gran- 
ular flows may be described through the dependency 
of the effective friction /i* (ratio of shear a to normal 
P stresses) and of the solid fraction v on the inertial 
number I = ^^/rn/P (a 2D equivalent of the deflnition 
given in Sec. [l|, where all the quantities are measured 
locally. The annular shear flows being heterogeneous, we 
measure the relations between the local quantities, i^{r), 
^*{r) = arg{r) I arr{r) and J(r) = 7(r) yjm/arr{r) (or 
i{r)/ y/orrir) in dimensionless unit). Each simulation 
provides dynamic dilatancy and friction laws in a range 
of inertial number. In the following, we try to analyze 
the granular material as a continuum, consequently we 
do not take into account the flve flrst layers where wall 
structuration effects are signiflcant (see Fig.[5]and Fig. [6]). 



1. Dynamic friction law 

In the inertial regime, for / > 0.02, ^* increases ap- 
proximately linearly with / and nearly independently of 
the geometry (Fig.jSK): 

f^* (1) - f^niin + bl , (5) 

with /ij„i„ ~ 0.26 and 6 ~ 1. The agreement with the 
dynamic friction law measured in the homogeneous plane 
shear geometry is excellent 0, 0, |4^. In contrast, for 
lower values of /, a deviation from this Hnear relation is 
observed, depending on the geometry (Fig. [Hb). The ef- 
fective friction becomes smaller than /imi„, and this devi- 
ation increases as Ri decreases, that is to say as the stress 
gradient increases. Reciprocally, as Ri increases, that is 
to say as the stress distribution becomes more homoge- 
neous, the results of the annular geometry tend to the 
ones of the plane shear geometry. This reveals that the 
simple relation between effective friction fj,* and inertial 
number / does not depend on the stress distribution in 
the inertial regime, and is then quite general (see 0, [s^l 
for flows down an inclined plane) , while it fails in the qua- 
sistatic regime. In plane shear, /ij^i„ may be considered 
as the internal friction in the critical state [H^l- This 
is the maximum value of fi* supported by the granular 
material, before it starts to flow quasistatically. With a 
heterogeneous stress distribution, the granular material 
is able to flow below this level. 

We call Xin the width of the inertial zone. Using 
Eqn. ^ and (O, we deduce that: 

X^n{Ve,Ri) - (^^S{Ve,R,)/l^;n^n~^^ R^- (6) 

We also conventionally deflne the width of the shear 
zone A/oc through vg{Ri + Xioc) = Vg/10. Fig. [9K and 
b shows Xin(Vg) and XiodVg) in geometry R50. We no- 
tice that Xin smoothly increases from zero with Vg , while 
Xioc seems to saturate at a low value for low Vg (global 
quasistatic regime) and at a high value for high Vg (this 
is related to an apparent velocity discontinuity near the 
wall, suggesting increasing collisional effects in the flrst 
layers), with a sudden increase for Vg between 0.3 and 1. 
We notice that in the experiment of [IBl, the shear zone 
invades the whole gap for high enough Vg. 

In a given geometry, for a small enough shear velocity 
Vg, the inertial zone disappears, and the whole system is 
in the quasistatic regime. Fig. [10] then shows again that 
the effective friction /i* is no more a function of /. 

2. Dynamic dilatancy law 

We observe a Hnear decrease of soHd fraction as a 
function of inertial number /, independently of the ge- 
ometry in the inertial regime (Fig. [11] and [T2|) . and v 
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FIG. 8: (Color online) Dynamic friction law (a) in linear 
scale (the solid line indicates a slope ~ 1^ for particle co- 
efficient of friction /i = and /i = 0.4. Dynamic friction 
law in semi-logarithmic scale for (h) = 0.4 and (c) fj, = 0. 
Different geometries: (U) i?25, (•) Rho,Jk) Rioo, (f) -R200. 
Ve = 2.5. Comparison with plane shear fji] (TJ). 



tends to a maximum value Vmax, which identifies to the 
solid fraction in the critical state. We can then write the 
dynamic dilatancy law: 

u{I) ~ Vrnax " 0,1, (7) 

with fmax — 0.82 and a ~ 0.37. The agreement with 
the dynamic dilatancy law measured in the homogeneous 
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FIG. 9: (a) Width of the inertial zone Xin as a function of 
Vg, as deduced from Eqn. ^ and Fig.\^. The solid line rep- 
resents the function: Xin = 50(^/1+051^^-1). (b) Width 
of the localization zone Xioc as a function of Ve. Geometry 
Rso- 
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FIG. 10: (Color online) Effective friction fj.* as a function of 
the mertial number I (ir) Ve = 0.00025, (j) Ve = 0.0025, 
(•)Ve ^ 0.025, (m)Ve^ 0.25, (*')Vb= 0.5, (^)Vb = 1.0, 
(■<*) Vg = 1.5, (a) Ve = 2.5. The solid line corresponds to 
/i* = 0.26. Geometry Rso- 



plane shear geometry is excellent 0, However, far 
from the walls, in the region where the material is less 
deformed and so remains in its initial dense state, higher 
values of v are observed. Fig. [12] also indicates that the 
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inner wall induces further dilation. 
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FIG. 11: (Color online) Dynamic dilatancy law (the solid line 
indicate a slope ~ —0.37^ for different geometries: (U) R25, 
(•) R50, (k) RiQo, (f)R200- Comparison with plane shear 
(D). 
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FIG. 12: (Color online) Dynamic dilatancy law in semi- 
logarithmic scale. Different geometries: (U) R25, (•) R50, 
(k) Rim, ('^ ) R200- Comparison with plane shear Jj/ (Oi). 
The dashed line corresponds to Umax = 0.82. 



3. Frictionless grains 

As shown in Fig. [8^ and Fig. [TTl the microscopic fric- 
tion coefficient ^ has a significant infiuence on the con- 
stitutive law parameters. Those figures also reveal good 
agreement with homogeneous shear simulations [1]. The 
solid fraction remains a linearly decreasing function of / 
(with a fast change in the quasistatic limit). The slope a 
is not affected, while Vmax increases to ~ 0.85. The dy- 
namic friction law keeps the same tendency but is shifted 
toward smaller values of friction. The linear approxima- 
tion with ~ 0.11 (Eqn. P) fails for / < 0.01. We 
notice that the range of validity of the dynamic friction 
law is much larger than for frictional grains, and that it 
does not seem to depend on the geometry. 



Those differences are likely related to some peculiar- 
ities of assemblies of frictionless grains The qua- 
sistatic limit, in such materials, is only approached for 
much smaller values of / than in the frictional case, and 
IJL^in is itself considerably lower. As a consequence on 
may expect a wider inertial zone. Moreover, as the criti- 
cal solid fraction coincides with the random close packing 
value no solid- like region of the system can be pre- 
vented from fiowing because of its density. 

4. Comparison with previous studies 

The validity of the constitutive law, once suitably gen- 
eralized to three dimensions, was successfully tested in 
fiows down a heap between lateral walls 0, In 
that case the velocity field, as deduced from numeri- 
cal computations in which the viscoplastic law was im- 
plemented, exhibits a more complex three-dimensional 
structure. Predicted velocities at the free surface agreed 
closely with experimental results. 

Thus, the applicability of the constitutive law as a re- 
lation between local values of non-uniform strain rate 
and stress fields, which we just established in 2D an- 
nular shear fiow, was previously checked in the 3D case 
of a laterally confined gravity-driven fiow. The validity 
of such an approach should be restricted to situations in 
which the characteristic length for stress or strain rate 
variations, say I, is significantly larger than the grain 
size. In annular shear, one has I = Ri, whereas the finite 
width w of the channel was found in 0, to control 
the gradients, I = w. As in units of grain diameters, 
varies between 25 and 200 here, while the interval of w 
extends between 16.5 and 500 in 0, similar levels of 
heterogeneity are explored. 

B. Internal variables 

We now discuss how internal variables, which profiles 
are discussed in App. C, scale with the inertial number /, 
revealing local state laws, consistent with the one mea- 
sured in homogeneous shear fiows. 

We observe a relation like Z = Zmax ~ el^ (with 
Zmax ~ 3) between coordination number Z and / on 
Fig. [in nearly independent of the geometry. 

We do not observe a general relation between the mo- 
bilization of friction M and /, but an asymptotic con- 
vergence for growing Ri toward a relation M w g/'' 
(Fig. [T4|) . For this quantity, there is no satisfactory agree- 
ment with the homogeneous shear case. 

We analyse the fiuctuations of orthoradial velocity Svg 
normalized by the natural scale 7 as a function of / 
(Fig. [T5| . In the quasistatic regime, the development 
of collective and intermittent motions (see [H, in an- 
nular shear and j6l| for a recent review) explain the 
increase of these relative fiuctuations. For higher values 
of the inertial number /, we observe that Sve/j 1. On 
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FIG. 13: (Color online) Coordination number Z as a function 
of the inertial number I (the solid line represents the function 
Z = 2.95 - T.eS/"'^^; for different geometries. (M) R25, (•) 
R50, (k) Rim, (f) -R200, (O) plane shear Q/. Vg = 2.5. 
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FIG. 15: (Color online) Relative fluctuations as a function of 
the inertial number I (the solid line represents the function 
Svg/i^d) = 1 + 0.077-0-^;. (M) R25, (•) R50, (A) Rioo, C^) 
R200, n plane shear Q/. Vg = 2.5. 
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where the shear stress at the inner wall S depends both 
on Vg and on Ri (see Sec. IIII Ap . As shown from the 
measurements drawn in Fig.fTBb. for a large value of Vg, 
S is high in small geometries and strongly decreases as Ri 
increases. We now integrate this relation over the range 
of validity of the dynamic friction law, this is to say in 
the inertial zone Ri Rin = Ri + Xin, from which we 
get: 



FIG. 14: (Color online) Mobilization of friction M as a func- 
tion of the inertial number I (the solid line represents the 
function M — 0.73/"'^^ for different geometries: (W) R25, (*) 
R50, (k) Rioo, (y) R200- Vg = 2.5. 



the whole, we propose to describe the dependency by the 
equation Svg/j = 1 + c/"'' (Fig. [15]). Experimental re- 
sults (2^. [26| show that Sv cc 7°-^. Dividing this relation 
by 7, we get an exponent equal to —0.6, close to exponent 
d= —0.7, deduced from the previous fit. 



vg{r) 



2br 



H — r m[r) + cr. 



(9) 



The constant c is determined by the value of the velocity 
at the inner wall, called V^ , which is smaller than Vg, 
reveahng some sliding at the wall as previously noticed 
(see Sec.iniB]). 



•-R^ \n{R,) + cR^ 



2b b 

On the whole, the velocity profile is equal to: 



(10) 



V. CONSEQUENCES FOR THE SHEAR 
LOCALIZATION AND THE MACROSCOPIC 
BEHAVIOR 

Using the constitutive law established in Sec. IIV| we 
now show that it is possible to understand some obser- 
vations described in Sec. IIIII 

Still using dimensionless units, since the pressure P is 
constant in the system and the shear stress is given by 
Eqn. ([3]), the dynamic friction law Eqn. ^ provides the 
following equation for the velocity profile vg{r): 



VB{r)^V+ — + 
Hi 




■In 



b "'\R, 

(11) 

An absolute measurement of Vg^ happens to be diffi- 
cult, considering the wall effect that disturbs the mate- 
rial behavior in a layer of a few grains near the inner wall 
(like in Fig. [5] for the shear rate 7). Consequently, we 
obtain this quantity from a fit, and a comparison with 
the measured velocity profiles is shown on Fig.fTBb. The 
agreement is excellent, suggesting once more the validity 



11 



of the dynamic friction law. The shding increases when 
Vg increases and, as shown in Fig. [T6b . increases when 
Ri decreases. 




FIG. 16: (Color online) (a) Influence of the geometry on (M) 
/Vo (fit) and (O) S (measurement) (Vg = 2.5). (b) Veloc- 
ity profiles: comparison between the measurements (U) R25, 
(») R30, (^) Rioo, (f) R200 and the prediction of Eqn. I|lip 

(solid lines). The velocity profiles are limited to the steady 
zone Ri ^ Ri + Rsteady (Ve — 2.5/ 



We now try to predict the S{Ve) relation, which was 
measured and fitted in Fig. [Sb. It is clear that in the 
global quasistatic Hmit, as Ve 0, S ^ l^min- We now 
write a boundary condition at the limit of the inertial 
zone Rin- Having used the dynamic friction law Eqn. |[5]), 
we necessarily have j{Rin) = 0, as appears in the fitted 
curve in Fig [T6b . Then, beyond Rin, if this dynamic 
friction law was still valid , 7(r) would be strictly equal to 
zero, so that vg{r) would be equal to Cr, with a constant 
C. The sole possibility is C = since the velocity must 
be equal to zero at the outer wall. We conclude that 
Vin = ve{Rin) ~ 0. This conclusion is wrong, as it is 
clear in Fig [TBb . and has already been discussed: the 
dynamic friction law fails in the quasistatic regime, and 
we shall come back to this point just after the discussion 
of the Siye) relation. The previous assumption writes: 



= 14^ 
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Since Eqn. ^ is equivalent to Rin/R-i = \/WJKi 
get the following implicit S{Vf^) relation: 



S 



1 - In 
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(13) 



For simplicity, we take V^ = Ve in the comparison with 
the measurements, drawn in Fig [17] for two geometries. 
The agreement is very satisfactory considering the previ- 
ous simplifying assumptions. With the increase of Ri, 
the difference between Vf^ and Vg decreases, explain- 
ing the better results for i?200- For small Vg, we write 
S{Vg) = ^^mini^ + fi^e))- A simple development gives 
/ ~ a/t-t^VW. For i?, = 50, / ~ 0.55VW, which is 
0.5V'g°-^'^ used in Fig. [3b. According 



close to the fit / 
to this analysis, S becomes proportional to Vg for much 
larger values, not usually accessible. 



0.60 



CO 




FIG. 17: (Color online) Shear stress at the wall S as a func- 
tion of the wall velocity Vg . Comparison between the measure- 
ments: (•) Rso and (J ) R200 and the predictions of Eqn. I|13p. 

The solid and dashed line respectively indicate the results for 
Rso and R200 . 



We now come back to the limit of the dynamic fric- 
tion law, Eqn. ((5|), in the quasistatic limit, as shown in 
Fig. [Hb- A large portion of the velocity profile in the 
steady quasistatic regime is shown in Fig. [T6b . As a first 
approximation, the velocity can be considered exponen- 
tial in this region, so that we write: 



vg{r) = Vin exp- 



r - R.. 



X 



qs 



(14) 



with Xqs the characteristic length in the quasistatic re- 
gion, measured in Fig. [TSj (which are slightly larger than 
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the one estimated for very small Vg, that is to say when 
the quasistatic zone invades the system), and Vin is not 
equal to zero contrarily to the previous simple approxi- 
mation but, using Eqn. ifTTj) to: 



2b 



S 



(15) 
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FIG. 18: (Color online) Characteristic length Xqs obtained 
from the velocity profiles in the quasistatic zone. Different 
geometries: (M) R25, (•) Rso, (k) Rum, (W ) i?200- Ve = 2.5. 



Still using dimensionless units, since the pressure P is 
constant in the system, we deduce that, for r > i?™, the 
inertial number is equal to: 



J(r) 



1 . n /r-i?, 

Vin CXp - 



(16) 



Since i?i„ ^ Ag,,, we may write: 



^(m ) - T— Gxp ^ 

■^qs -^qs 



from which we obtain: 
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(17) 



(18) 



This prediction is in close agreement with the measure- 
ments, as shown in Fig. [191 



VI. CONCLUSION 

We first summarize the results presented in this paper, 
before discussing the questions raised by those conclu- 
sions. 
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FIG. 19: Effective friction n* as function of the inertial num- 
ber I , in the quasistatic zone. Comparison between the mea- 
surements: (M) R25, (•) R50, (k) Rim, ('^ ) R200, (O) plane 
shear 07 and the prediction of Eqn. (1 18 II (solid lines). The 
dashed line indicates ^I'min ~ 0.26. Vg — 2.5. 



As described in Sec. m we have studied through dis- 
crete simulations steady annular shear flows of a model 
granular material, made of a slightly polydisperse assem- 
bly of frictional dissipative disks, prescribing the rotation 
rate of the inner wall and the pressure exerted by the 
outer wall, and varying dimensionless shear velocity Ve 
and size Ri of the system. 

The first step (Sec. lIIip has consisted in measuring var- 
ious quantities, either global as the dimensionless shear 
stress at the inner wall as a function of Ve, or local 
as the profiles of stress components, velocity, solid frac- 
tion and some internal variables (coordination number, 
mobilization of friction, velocity fluctuations, shown in 
App. C). This has allowed to distinguish, at the global 
scale, that is to say as a function of Ve, between rate 
dependent and rate independent behaviors. 

Inspired by our previous rheophysical analysis of ho- 
mogenous shear flows of disks the second step 
(Sec. IIVP has explored the validity of constitutive law 
for inertial regime if applied locally in such an hetero- 
geneously sheared material. We have shown that the 
dynamic friction and dilatancy laws observed in homo- 
geneous shear flows are exactly recovered, when using 
the local state parameter / called inertial number. Scal- 
ing laws for internal variables as function of / have also 
been observed. This analysis has clearly distinguished an 
inertial zone close to the inner wall where the constitu- 
tive law is relevant and a quasistatic zone away from it, 
where it fails. 

The last step (Sec. |V| has explained how it is possible 
to predict some observations presented in the first step, 
when using the inertial constitutive law identified in the 
second step. We have focused on two basic quantities, 
which are most often discussed in the studies of annular 
shear fiows of granular materials, the macroscopic S{Ve) 
relation and the microscopic velocity profiles. The satis- 
factory agreement between the prediction and the mea- 



13 



surements should not be surprising, considering the sec- 
ond step. However this analysis has precisely pointed out 
two important issues, related to boundary conditions in 
such a heterogeneously sheared system, one at the shear- 
ing wall, and the other at the transition between inertial 
and quasistatic zone. 

Close to the inner wall, as previously shown in different 
configurations, various quantities (solid fraction, ratio of 
normal stresses, rotation velocity...) present singular be- 
haviors. The translation velocity reveals significant slid- 
ing for sufficiently high Ve , even with the large roughness 
used in this study. We have shown that the value Vg^ of 
this sliding velocity is an important ingredient for a good 
prediction. This means that a detailed understanding of 
the rheophysics of the granular materials in the very first 
layers near a rough wall is of great relevance. Apart from 
the characteristics of the granular material itself, the rel- 
ative infiuences of Ve, Ri and the wall roughness must be 
taken into account. Comparisons between physical ex- 
periments and discrete simulations are described in [43 |. 
Considering the frustration of the particle rotation im- 
posed by the wall, Cosserat models might be adapted 
to describe this interface zone [g^I, as done by (iol . [63| 
for annular shear. Another discussion of the boundary 
condition at the wall is proposed in (63 |. 

The transition between inertial and quasistatic zones 
is a second puzzHng issue. Considering the constitutive 
law identified in homogeneous shear fiows, the granular 
material should reach the so-called critical state in the 
quasistatic limit (when / — > 0), in which it fiows rate 
independently with an effective friction tan (f) and a soHd 
fraction Vc- Beyond this limit (for S/P < tan0 and/or 
V > Vc), the granular material, being in a solid-like state, 
should not be able to fiow. However (apparently un- 
bounded) creep fiows are observed in this nominally soHd 
regime. This creeping behavior is well known in free sur- 
face fiows, where an exponential velocity profile has been 
clearly evidenced with a characteristic length of the or- 
der of one grain diameter [sl, [13, [13, [l^l . In the annular 
shear geometry, a similar behavior is observed but the 
characteristic length increases as Ri increases, that is to 
say as the stress gradient decreases, or as the stress field 
becomes more homogeneous. 

For a sufficiently small Ve, there is no more inertial 
zone, so that both boundary conditions occur at the same 
place, the inner wall. Considering the typical values of 
the parameters in the systems which have been studied 
experimentally or through discrete simulations, we notice 
that this corresponds to the usual case. The understand- 
ing of such a situation merges the two previous problems: 
the behavior of a granular material close to an interface 
and in the quasistatic regime, together with the hetero- 
geneity of the stress field. 

The already noticed observation of collective and in- 
termittent motions in this quasistatic regime has driven 
the development of several rheological models (see jH, 
[sol . [gqI . [toI and [6l| for a recent review) : diffusion equa- 
tion for the fiuctuations, transmission of forces at the 



scale of correlated clusters, two-phase fiuid model with 
order parameter, activation of rearrangements through 
the fiuctuations of velocity or forces, occurring either at 
the boundary of the inertial zone, or at the inner wall in 
the global quasistatic limit. 

Our understanding is far from complete and requires 
further studies, merging physical experiments, discrete 
simulations and theoretical developments. For instance, 
we have not measured the fabric in the quasistatic zone, 
although its importance has been clearly evidenced in 
homogeneous shear [H^. [53|. We have not discussed the 
evolution of the internal variables in the transient regime 
(evolution from initial to steady state), or in a shear re- 
versal regime (25l . [3l| , as should be qualitatively possible 
using simplified microscopic description [7l| . We have re- 
stricted our attention to velocity controlled shear fiows, 
so that it was not possible to study the fiow threshold. A 
specific study of the jamming mechanisms should be per- 
formed by controlling the shear stress jl^l- We have not 
discussed the infiuence of the r oug hness on the interface 
behavior, for which we refer to [42|. We may also wonder 
to what extent the conclusions drawn for granular mate- 
rials differ for other complex fiuids made of interacting 
elements (dense suspensions, foam, emulsions...) (7^.[7^. 



Appendix A : Periodic boundary condition 

Each grain the center of which is in (r, 6*) with < 9 < 
& is associated to a collection of copies with centers in 
r,9 + kO where k is an integer. The corresponding ve- 
locities, accelerations and forces are related by rotations 
of angles kQ. 

Every time a grain moves out of the simulation cell, 
one of its copies moves in by the opposite boundary, sim- 
ilarly to the usual case of periodic boundary conditions 
by translation. However the velocities, accelerations and 
forces are affected by a rotation of ±9. 

The situation of the contact of two grains i and j where 
9i is close to Q and 9j is close to zero is described in 
Fig. EOl More precisely, i is in contact with the copy 
j' of j, obtained by rotation of an angle Q, while j is 
in contact with i' obtained by rotation of i of an angle 
—9. To evaluate the forces acting over grain i we have to 
use the normal and tangential unit vectors n^' (pointing 
from i to j') and ty' (such that {nij' ,tij') is positively 
oriented), respectively, and the motion of the grain j' , 
while for j we have to use corresponding fljir and tjii 
and the motion of i'. Vector n^' is not, as usually, equal 
to — fiji' , but to its image obtained by rotation of an angle 
-9. 

We have measured the infiuence of the periodic bound- 
ary condition comparing the radial profiles of various 
quantities as a function of 9 (tt/IG, tt/S, 7r/4, 7r/2, tt 
and the whole ring 27r) for the geometry Ri = 25 and 
i?o = 50. As an example, we show on Fig.[2T]the profiles 
of the orthoradial velocity. As expected, the results are 
all the more consistent as the value of 9 increases. In 
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FIG. 20: Periodic boundary conditions. 



this case, 6 ~ 7r/2 already gives a very good result. 




FIG. 21: ( Color online) Velocity profiles vg{r) /Ve for different 
values o/e (rad). (•)& = n/lQ, (k) Q = tt/S, ) Q = n/2, 
(m)@ = 2-K. Ri^ 25, Ro = 50, Vg = 2.5. 



We quantify the deviations of the velocity profiles V0{r) 
by means of an indicator of relative error. The velocity 
tends to zero as the distance from the inner wall. To avoid 
inconsistencies due to values close to zero in the frame 
of the usual definition of relative error, and to give more 
weight to the values close to the inner wall, we propose 
to calculate the relative error over variable Fq (r) = — 
■ye(r, 9) : 



£(6) = 



Ro — Ri 



iRi 



, Feir) - F2eir) 
' F2e(r) 



dr. 



(19) 



e(9) is simply the sum over the whole geometry of the 
relative error of the variable F for a certain value of Q 
compared to the result for a system twice as large (29). 

On Fig. [22k . we observe a clear decrease of the error 
indicator e(9) as we increase the value of 9 for the small- 
est geometry {Ri = 25 and Ro = 50). The same analysis 



for a larger geometry {Ri = 100 and Ro = 200) shows 
better results for smaller values of 9. This shows that 
the influence of 9 on the results depends on the size of 
the system. We try to relate both parameters in Fig. [22b. 
where we plot e(9) as a function of the angular sector 
length at the inner wall {SRi). We observe that a good 
accuracy of the results can be achieved with a length 
QRi > 40 for geometries with Ri > 25. Based on this 
consideration, we have chosen the values of 9 for each of 
our geometries (Tab. I). 
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FIG. 22: (Color online) Relative error e on the orthoradial 
velocity vg (Vg — 2.5), (a) as function ofQ, (h) as function 
of the inner wall length QRi. (W) Ri = 2b and Ro ~ 50, (•) 
Ri = 100 and Ro = 200. 



Appendix B : Averaging method 

Considering the revolution symmetry of our system, 
the radial profiles of different quantities (orthoradial ve- 
locity ve{r), coordination number Z(r), etc.) are ob- 
tained by an averaging procedure over coordinate along 
the coordinate r (Fig. [23|) . To each of the n grains i are 
associated different scalar quantities G". We define a 
weight function i/'i(r) as the intercept angle defined on 
Fig.|2a(cos(V'«(r)/2) = {r"^ +rf - df/ 4.) / {2rr{) for a disk 
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of diameter di). Some variables, like solid fraction v, are 
averaged over the whole space, while others, like the co- 
ordination number Z, have no sense outside the grain 
space. This leads to the two following definitions of the 



average : 



and : 



1 " 



(20) 



The stress tensor of each grain is defined according 
to [lil (with A"" = TTdj /4 the grain area) : 



(24) 



The first term is associated to the contact forces, and 
the second one to the velocity fiuctuations. The radial 
profiles of the components of the stress tensor are : 



< G >' (r) = 



(21) 




FIG. 23: Various quantities associated to a grain i. 



Applying this principle, we determine the solid fraction 
profile fir) as follows 



1 

i— 1 



(22) 



where the value of f ^ is naturally equal to 1 . This means 
that < G > and < G >' are simply related by the solid 
fraction: < G >=< ly >< G >'. 

We take into account the variation of vectorial and 
tensorial quantities inside the grains, when written in 

the polar basis e^) = ) and ee = 'J^'^^ 

Hence, the radial profiles of the velocity components are 



n 

E 

1=1 



(23) 



^"/3(o = E 



.{<!>)■ a! ■ ep{cl,)dcl,. (25) 



Since we try to analyze the granular material as a 
continuum (except for the very first layers near the wall) , 
we consider the coarse-grained variations of the quan- 
tities by smoothing the profiles through central moving 
averages of 3d length (if not otherwise indicated) . The re- 
maining fiuctuations would disappear with an increase of 
the simulation time At over which the data are averaged. 



Appendix C : Internal variables 

Coordination number Z is the average number of con- 
tacts per grain. In the inertial regime, the general ten- 
dency is a decrease of Z as the shear rate 7 increases 
(that is to say for increasing Vg in Fig. [21]) . For smaller 
values of 7 (corresponding to smaller values of Ve or to 
a larger distance from the inner wall) the coordination 
number Z approaches a limiting value, slightly above 3. 
Such a limit is in rough agreement with other numeri- 
cal observations of the critical state of frictional disks. 
Ref. [H^l thus reports Z ~ 3.6. The somewhat lower val- 
ues observed in our case are likely due to the larger strain 
rates, and to the remaining infiuence, on the quasistatic 
region of limited width, of the more agitated inner zone. 

We define the mobilization of friction as ratio M = 
Zs/Z, where Zs is the average number of sliding contacts 
per grain 0, [73|. Fig. [25] shows that M increases as the 
shear rate increases, whether through an increase of Ve or 
a decrease of the distance from the inner wall. We notice 
that the stabilization of the M{r) profile occurs for Ve < 
0.0025, a value much smaller than the one required for the 
stabilization of the other studied quantities {Ve < 0.025). 

For any quantity q{r) averaged in space (along 9) and 
in time, we may define its fiuctuation: 



qir,e)^d9-qirf, 



(26) 



where q{r, 6*)^ is averaged in time. We measure the 
fiuctuations of the translational and rotational velocities 
Sve{r), 5vr{r) and 5w(r). Our analysis (long time scale) 



16 




r-R. 

I 

PIG. 24: (Color online) Influence of shear velocity Vg on the 
coordination number profiles Z{r). (J)V0 = 0.0025, (•)Vg — 
0.025, (m) Ve = 0.25, (t^) Ve = 0.5, (^) Vg = 1.0, (a) 
Vg = 1.5, (k)Vg = 2.5. Geometry Rso- 



10° 




r-R. 

I 

FIG. 25: (Color online) Influence of shear velocity Vg on the 
mobilization of friction profiles M{r) . (ic ) Vg = 0, 00025, (W ) 
Vg = 0.0025, (•)Vg = 0.025, (M) Vg = 0.25, (k) Vg = 2.5. 
Geometry Rso- 



takes into account both the small fluctuations around the 
mean motion (in the cage formed by the nearest neigh- 
bors), and the large fluctuations associated to collective 
motions [t^ . 

Fig. [26] flrst shows that the general amplitude of the 
fluctuations increases with Vg. Then, for various Vg, they 
reveal a strong decay of the fluctuating quantities close 
to the inner wall, comparable to that of the respective 
average quantities, consistently with previous observa- 
tions [1, [l3, [13, S, m, [Hi. This decay is still true at 
larger distances for Svg and Su> (with an increase close to 
the outer wall). We also notice a stabilization of Svr, 
which occurs at r — Ri sa 10 for Vg = 0.025 and at 
r — i?i w 20 for Vg — 2.5, that is to say precisely when the 
solid fraction v reaches a value w 0.82 (Fig. [6]). Above 
this critical value of i^, the material would be so com- 



pact that the radial motions would take place as a block. 
Fig. [26k shows the equality of Svg and Svr before the sta- 
biHzation of 5vr{r), while Fig. [26b shows the systematic 
equality of Svg and Suj/2. 




r-R. 

I 

FIG. 26: (Color online) (a) Comparison between the profiles 
of the fluctuations of the radial velocity Svr (hollow symbols) 
and of the orthoradial velocity Svg (full symbols), (b) Com- 
parison between the profiles of the fluctuations of the angular 
velocity 5lo/2 (hollow symbols) and of the orthoradial veloc- 
ity Svg (full symbols) for different wall velocities Vg: ('V, W) 
Vg = 0.0025, (o,») Vg = 0.025, (D.M) Vg = 0.25, (A, A) 
Vg = 2.5. Geometry R50. 
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